load the libraries
library("readxl")
library("tidyverse")
library("rstatix")
library("DESeq2")
library("EnhancedVolcano")
library("ggplot2")
library("pmp")
library("ggpubr")
library("ComplexHeatmap")
library("limma")
library("qgraph")
load the files and assign them to respective variable names.
lipid <-as.data.frame(read_excel("italian cohort lipid study.xlsx",sheet = 2))
lipid<-lipid[-1,-277]
clinical_data<-lipid[,1:4]
lipid_data<-lipid[,5:ncol(lipid)]
annotations <- as.data.frame(read_excel("annotations.xlsx",col_names =T))
colnames(lipid_data)<-make.unique(annotations$lipid)
rownames(lipid_data)<-clinical_data$sample_ID
lipid_data<-as.data.frame(lipid_data)
summary(lipid_data[,1:10])
## CE_14:0 CE_15:0 CE_16:0 CE_16:0.1
## Min. :0.0000000 Min. :1.885e-05 Min. :0.0001391 Min. :0.003984
## 1st Qu.:0.0003600 1st Qu.:8.394e-05 1st Qu.:0.0002549 1st Qu.:0.006530
## Median :0.0004367 Median :1.192e-04 Median :0.0004512 Median :0.007532
## Mean :0.0004403 Mean :1.187e-04 Mean :0.0005237 Mean :0.007500
## 3rd Qu.:0.0005176 3rd Qu.:1.476e-04 3rd Qu.:0.0006468 3rd Qu.:0.008625
## Max. :0.0008344 Max. :2.630e-04 Max. :0.0016558 Max. :0.010846
## CE_16:1 CE_16:1.1 CE_17:0 CE_17:1
## Min. :0.000e+00 Min. :0.001046 Min. :2.787e-05 Min. :0.0000000
## 1st Qu.:3.230e-05 1st Qu.:0.002081 1st Qu.:7.312e-05 1st Qu.:0.0001557
## Median :5.871e-05 Median :0.002471 Median :1.012e-04 Median :0.0002392
## Mean :8.938e-05 Mean :0.002615 Mean :1.105e-04 Mean :0.0002002
## 3rd Qu.:1.196e-04 3rd Qu.:0.003015 3rd Qu.:1.453e-04 3rd Qu.:0.0002709
## Max. :4.621e-04 Max. :0.006832 Max. :2.295e-04 Max. :0.0003762
## CE_18:1 CE_18:1.1
## Min. :0.0003304 Min. :0.01170
## 1st Qu.:0.0005593 1st Qu.:0.02253
## Median :0.0008584 Median :0.02588
## Mean :0.0009193 Mean :0.02585
## 3rd Qu.:0.0011769 3rd Qu.:0.02856
## Max. :0.0023625 Max. :0.04234
table(clinical_data$NAFLD)
##
## 0 1
## 120 21
scale the data
lipid_scaled<-scale(lipid_data)
summary(lipid_scaled[,1:10])
## CE_14:0 CE_15:0 CE_16:0 CE_16:0.1
## Min. :-3.11746 Min. :-2.21915 Min. :-1.0770 Min. :-2.35097
## 1st Qu.:-0.56875 1st Qu.:-0.77224 1st Qu.:-0.7526 1st Qu.:-0.64815
## Median :-0.02568 Median : 0.01217 Median :-0.2030 Median : 0.02164
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.54740 3rd Qu.: 0.64303 3rd Qu.: 0.3445 3rd Qu.: 0.75223
## Max. : 2.79007 Max. : 3.20715 Max. : 3.1696 Max. : 2.23800
## CE_16:1 CE_16:1.1 CE_17:0 CE_17:1
## Min. :-1.0552 Min. :-1.9408 Min. :-1.7820 Min. :-1.8814
## 1st Qu.:-0.6738 1st Qu.:-0.6612 1st Qu.:-0.8060 1st Qu.:-0.4180
## Median :-0.3621 Median :-0.1781 Median :-0.2014 Median : 0.3673
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.3566 3rd Qu.: 0.4948 3rd Qu.: 0.7513 3rd Qu.: 0.6645
## Max. : 4.3995 Max. : 5.2148 Max. : 2.5666 Max. : 1.6547
## CE_18:1 CE_18:1.1
## Min. :-1.3858 Min. :-2.662815
## 1st Qu.:-0.8471 1st Qu.:-0.625705
## Median :-0.1432 Median : 0.004171
## Mean : 0.0000 Mean : 0.000000
## 3rd Qu.: 0.6064 3rd Qu.: 0.509735
## Max. : 3.3961 Max. : 3.101539
conduct PCA analysis
res_pca_1<-prcomp(lipid_scaled,centre=T,scale=F)
res_pca <- as.data.frame(res_pca_1$x)
res_pca$steatosis<-clinical_data$NAFLD
res_pca$sample<-clinical_data$sample_ID
theme<-theme(panel.background = element_blank(),panel.border=element_rect(fill=NA),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),strip.background=element_blank(),axis.text.x=element_text(colour="black"),axis.text.y=element_text(colour="black"),axis.ticks=element_line(colour="black"),plot.margin=unit(c(1,1,1,1),"line"))
percentage <- round(res_pca_1$sdev / sum(res_pca_1$sdev) * 100, 2)
percentage <- paste( colnames(res_pca), "(", paste( as.character(percentage), "%", ")", sep="") )
plot1<-ggplot(res_pca,aes(x=PC1,y=PC2,label= sample,color=steatosis ))
plot1<-plot1+geom_point(size=4)+theme+ xlab(percentage[1]) + ylab(percentage[2])
print(plot1)
rename the sample in clinical data
clinical_data$NAFLD[clinical_data$NAFLD==0]<-"steatosis0"
clinical_data$NAFLD[clinical_data$NAFLD==1]<-"steatosis1"
head(clinical_data)
## sample_ID patient_no NAFLD sex
## 2 N1 1 steatosis0 F
## 3 N2 1 steatosis0 F
## 4 N3 1 steatosis0 F
## 5 N4 1 steatosis0 F
## 6 N5 1 steatosis0 F
## 7 N6 1 steatosis0 F
Extract the steatosis 0 and steatosis 1 sample separately
clinical_data_steatosis0<-subset(clinical_data, NAFLD=="steatosis0",select=sample_ID:sex)
clinical_data_steatosis1<-subset(clinical_data, NAFLD=="steatosis1",select=sample_ID:sex)
head(clinical_data_steatosis1)
## sample_ID patient_no NAFLD sex
## 76 N75 15 steatosis1 F
## 77 N76 15 steatosis1 F
## 82 N81 18 steatosis1 F
## 83 N82 18 steatosis1 F
## 86 N85 20 steatosis1 M
## 87 N86 20 steatosis1 M
head(clinical_data_steatosis0)
## sample_ID patient_no NAFLD sex
## 2 N1 1 steatosis0 F
## 3 N2 1 steatosis0 F
## 4 N3 1 steatosis0 F
## 5 N4 1 steatosis0 F
## 6 N5 1 steatosis0 F
## 7 N6 1 steatosis0 F
lets see how many steatosis0 and steatosis1 samples are present
nrow(clinical_data_steatosis0)
## [1] 120
nrow(clinical_data_steatosis1)
## [1] 21
There are more steatosis0 than steatosis 1 , hence lets conduct differential expression of lipid analysis between these two groups by subseting the samples into separate batches.
create 6 batches of distinct (without sample replacement of steatosis0) and conduct differential expression of lipid for each batch of steatosis0 and steatosis1 such as (20 vs 21) and identify the repeatedly arising significant gene between them.
Lets create 6 batches
rownames(clinical_data_steatosis0)<-1:120
batch1_clinical_data_steatosis0<-clinical_data_steatosis0[1:20,]
batch2_clinical_data_steatosis0<-clinical_data_steatosis0[21:40,]
batch3_clinical_data_steatosis0<-clinical_data_steatosis0[41:60,]
batch4_clinical_data_steatosis0<-clinical_data_steatosis0[61:80,]
batch5_clinical_data_steatosis0<-clinical_data_steatosis0[81:100,]
batch6_clinical_data_steatosis0<-clinical_data_steatosis0[101:120,]
now create the clinical data for those 6 batches by merging clinical batch of steatosis0 and steatosis1
batch1_clinical_data<-rbind(batch1_clinical_data_steatosis0,clinical_data_steatosis1)
batch2_clinical_data<-rbind(batch2_clinical_data_steatosis0,clinical_data_steatosis1)
batch3_clinical_data<-rbind(batch3_clinical_data_steatosis0,clinical_data_steatosis1)
batch4_clinical_data<-rbind(batch4_clinical_data_steatosis0,clinical_data_steatosis1)
batch5_clinical_data<-rbind(batch5_clinical_data_steatosis0,clinical_data_steatosis1)
batch6_clinical_data<-rbind(batch6_clinical_data_steatosis0,clinical_data_steatosis1)
extract the expression data for each of the batches now
batch1_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch1_clinical_data$sample_ID,]
batch2_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch2_clinical_data$sample_ID,]
batch3_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch3_clinical_data$sample_ID,]
batch4_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch4_clinical_data$sample_ID,]
batch5_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch5_clinical_data$sample_ID,]
batch6_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch6_clinical_data$sample_ID,]
conduct differential lipid expression analysis for batch-1
design_batch1 <- model.matrix(~0+batch1_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch1) <- c("Steatosis0","Steatosis1")
lipid_data_limma1<-t(batch1_lipid_scaled)
fit_data1 <- lmFit(lipid_data_limma1, design_batch1)
contrasts1 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch1)
fit2_data1 <- contrasts.fit(fit_data1, contrasts1)
fit2_data1<- eBayes(fit2_data1)
full_results1 <- topTable(fit2_data1, number=Inf,adjust = "BH")
full_results1 <- cbind(symbols = rownames(full_results1), full_results1)
rownames(full_results1) <- 1:nrow(full_results1)
lets view the plot for batch 1 and extract the significant genes of them
p_cutoff <- 0.05
fc_cutoff <- 1
topN <-10
full_results1 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch1<-full_results1[full_results1$adj.P.Val<p_cutoff,]
head(significant_batch1)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 TG_53:3 1.781418 0.007213676 10.84643 1.509766e-14 3.931815e-12 22.91345
## 2 DG_34:3 2.125923 -0.021706774 10.56533 3.714068e-14 3.931815e-12 22.01915
## 3 DG_43:6 -1.832345 0.029401360 -10.51652 4.346753e-14 3.931815e-12 21.86286
## 4 TG_56:4 2.004968 0.186773092 10.33753 7.759063e-14 3.931815e-12 21.28707
## 5 TG_56:5 2.036917 0.106174432 10.33529 7.815784e-14 3.931815e-12 21.27983
## 6 TG_50:3 1.766674 -0.181929281 10.30776 8.547425e-14 3.931815e-12 21.19090
dim(significant_batch1)
## [1] 191 7
conduct differential lipid expression analysis for batch-2
design_batch2 <- model.matrix(~0+batch2_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch2) <- c("Steatosis0","Steatosis1")
lipid_data_limma2<-t(batch2_lipid_scaled)
fit_data2 <- lmFit(lipid_data_limma2, design_batch2)
contrasts2 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch2)
fit2_data2 <- contrasts.fit(fit_data2, contrasts2)
fit2_data2<- eBayes(fit2_data2)
full_results2 <- topTable(fit2_data2, number=Inf,adjust = "BH")
full_results2 <- cbind(symbols = rownames(full_results2), full_results2)
rownames(full_results2) <- 1:nrow(full_results2)
volcano plot for batch 2
full_results2 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch2<-full_results2[full_results2$adj.P.Val<p_cutoff,]
head(significant_batch2)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 DG_38:5 -1.833957 0.008044289 -11.115877 3.355189e-15 9.260321e-13 24.37448
## 2 DG_43:6 -1.834214 0.030313308 -10.684875 1.383043e-14 1.908599e-12 22.96750
## 3 PI_39:4 -1.693211 0.032347451 -10.212007 6.717783e-14 6.180360e-12 21.39697
## 4 PS-O_38 -1.738166 0.021359607 -10.010747 1.327170e-13 9.157471e-12 20.72026
## 5 TG_56:4 1.928321 0.224161721 9.829341 2.461660e-13 1.358836e-11 20.10622
## 6 TG_53:3 1.766265 0.014605727 9.568558 6.022766e-13 2.770473e-11 19.21690
dim(significant_batch2)
## [1] 189 7
conduct differential lipid expression for batch 3
design_batch3 <- model.matrix(~0+batch3_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch3) <- c("Steatosis0","Steatosis1")
lipid_data_limma3<-t(batch3_lipid_scaled)
fit_data3 <- lmFit(lipid_data_limma3, design_batch3)
contrasts3 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch3)
fit2_data3 <- contrasts.fit(fit_data3, contrasts3)
fit2_data3<- eBayes(fit2_data3)
full_results3 <- topTable(fit2_data3, number=Inf,adjust = "BH")
full_results3 <- cbind(symbols = rownames(full_results3), full_results3)
rownames(full_results3) <- 1:nrow(full_results3)
volcano plot for batch 3
full_results3 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch3<-full_results2[full_results3$adj.P.Val<p_cutoff,]
head(significant_batch3)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 DG_38:5 -1.833957 0.008044289 -11.115877 3.355189e-15 9.260321e-13 24.37448
## 2 DG_43:6 -1.834214 0.030313308 -10.684875 1.383043e-14 1.908599e-12 22.96750
## 3 PI_39:4 -1.693211 0.032347451 -10.212007 6.717783e-14 6.180360e-12 21.39697
## 4 PS-O_38 -1.738166 0.021359607 -10.010747 1.327170e-13 9.157471e-12 20.72026
## 5 TG_56:4 1.928321 0.224161721 9.829341 2.461660e-13 1.358836e-11 20.10622
## 6 TG_53:3 1.766265 0.014605727 9.568558 6.022766e-13 2.770473e-11 19.21690
dim(significant_batch3)
## [1] 171 7
conduct differential lipid expression for batch 4.
design_batch4 <- model.matrix(~0+batch4_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch4) <- c("Steatosis0","Steatosis1")
lipid_data_limma4<-t(batch4_lipid_scaled)
fit_data4 <- lmFit(lipid_data_limma4, design_batch4)
contrasts4 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch4)
fit2_data4 <- contrasts.fit(fit_data4, contrasts4)
fit2_data4<- eBayes(fit2_data4)
full_results4 <- topTable(fit2_data4, number=Inf,adjust = "BH")
full_results4 <- cbind(symbols = rownames(full_results4), full_results4)
rownames(full_results4) <- 1:nrow(full_results4)
volcano plot for batch 4
full_results4 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch4<-full_results4[full_results4$adj.P.Val<p_cutoff,]
head(significant_batch4)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 DG_38:5 -1.583231 -0.22742195 -7.799571 4.108759e-10 5.608425e-08 12.923874
## 2 DG_34:3 1.536189 0.47264808 7.771243 4.539084e-10 5.608425e-08 12.826047
## 3 DG_43:6 -1.569167 -0.22818895 -7.687448 6.096114e-10 5.608425e-08 12.536394
## 4 PS-O_38 -1.269659 -0.29150250 -7.349281 2.012102e-09 1.388351e-07 11.363847
## 5 Oxochol 1.285460 -0.05608858 6.472756 4.516807e-08 2.493277e-06 8.312214
## 6 PI_41:6 -1.219338 -0.22827882 -6.095570 1.720613e-07 7.914820e-06 7.003015
dim(significant_batch4)
## [1] 129 7
Conduct differential lipid expression for batch 5
design_batch5 <- model.matrix(~0+batch5_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch5) <- c("Steatosis0","Steatosis1")
lipid_data_limma5<-t(batch5_lipid_scaled)
fit_data5 <- lmFit(lipid_data_limma5, design_batch5)
contrasts5 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch5)
fit2_data5 <- contrasts.fit(fit_data5, contrasts5)
fit2_data5<- eBayes(fit2_data5)
full_results5 <- topTable(fit2_data5, number=Inf,adjust = "BH")
full_results5 <- cbind(symbols = rownames(full_results5), full_results5)
rownames(full_results5) <- 1:nrow(full_results5)
volcano plot for batch 5
full_results5 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch5<-full_results5[full_results5$adj.P.Val<p_cutoff,]
head(significant_batch5)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 DG_36:2.1 -1.1993450 0.49101262 -5.958626 4.239642e-07 6.882647e-05 6.304671
## 2 PC_36:4 1.5397202 0.08691141 5.910031 4.987425e-07 6.882647e-05 6.147818
## 3 TG_53:2 -0.9889126 0.87228183 -4.635981 3.318158e-05 2.184353e-03 2.109665
## 4 PC_38:4 1.1610951 -0.12761206 4.623782 3.450870e-05 2.184353e-03 2.072151
## 5 TG_53:3 -0.7102970 0.86506896 -4.510687 4.957409e-05 2.184353e-03 1.725896
## 6 TG_53:6 -1.2651848 0.69742118 -4.465661 5.722572e-05 2.184353e-03 1.588845
dim(significant_batch5)
## [1] 42 7
Conduct different lipid expression for batch 6
design_batch6 <- model.matrix(~0+batch6_clinical_data$NAFLD,ref="steatosis0")
colnames(design_batch6) <- c("Steatosis0","Steatosis1")
lipid_data_limma6<-t(batch6_lipid_scaled)
fit_data6 <- lmFit(lipid_data_limma6, design_batch6)
contrasts6 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch6)
fit2_data6 <- contrasts.fit(fit_data6, contrasts6)
fit2_data6 <- eBayes(fit2_data6)
full_results6 <- topTable(fit2_data6, number=Inf,adjust = "BH")
full_results6 <- cbind(symbols = rownames(full_results6), full_results6)
rownames(full_results6) <- 1:nrow(full_results6)
volcano plot for batch 6
full_results6 %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_batch6<-full_results6[full_results6$adj.P.Val<p_cutoff,]
head(significant_batch6)
## symbols logFC AveExpr t P.Value adj.P.Val B
## 1 PC-O_44 1.4429320 0.10057341 5.973130 4.210383e-07 0.0001069434 6.326011
## 2 TG_54:2 -0.9965229 0.21923572 -5.685326 1.095100e-06 0.0001069434 5.405205
## 3 PC-O_36 1.3393702 -0.06098412 5.667329 1.162428e-06 0.0001069434 5.347761
## 4 PS-P_28 0.9521793 -0.14752453 5.376447 3.041048e-06 0.0002098323 4.422556
## 5 PC-O_40 1.4239131 0.06407383 5.038354 9.218303e-06 0.0005088503 3.357679
## 6 DG_36:2 -1.0665274 0.37025704 -4.825684 1.839361e-05 0.0008461062 2.695810
dim(significant_batch6)
## [1] 57 7
Find the intersecting significant lipids among the 6 batches
lipids_batch1<-significant_batch1$symbols
lipids_batch2<-significant_batch2$symbols
lipids_batch3<-significant_batch3$symbols
lipids_batch4<-significant_batch4$symbols
lipids_batch5<-significant_batch5$symbols
lipids_batch6<-significant_batch6$symbols
intersected_lipids<-Reduce(intersect,list(lipids_batch1,lipids_batch2,lipids_batch3,lipids_batch4,lipids_batch5,lipids_batch6))
print(intersected_lipids)
## [1] "TG_53:3" "DG_34:1" "TG_52:3" "TG_52:2.1" "TG_53:2"
## [6] "TG_53:7" "TG_54:2" "DG_36:2" "TG_53:6" "CE_16:0.1"
## [11] "Cholesterol"
conduct different expression of lipid only with the intersected lipids
common_lipid_scaled<-lipid_scaled[,colnames(lipid_scaled) %in% intersected_lipids]
design_final <- model.matrix(~0+clinical_data$NAFLD,ref="steatosis0")
colnames(design_final) <- c("Steatosis1","Steatosis0")
lipid_data_limma_final<-t(common_lipid_scaled)
fit_data_final <- lmFit(lipid_data_limma_final, design_final)
contrasts_final <- makeContrasts(Steatosis1 - Steatosis0, levels=design_final)
fit2_data_final <- contrasts.fit(fit_data_final, contrasts_final)
fit2_data_final<- eBayes(fit2_data_final)
full_results_final <- topTable(fit2_data_final, number=Inf,adjust = "BH")
full_results_final <- cbind(symbols = rownames(full_results_final), full_results_final)
rownames(full_results_final) <- 1:nrow(full_results_final)
print(full_results_final)
## symbols logFC AveExpr t P.Value adj.P.Val
## 1 TG_53:3 -1.0295330 3.430077e-17 -4.531666 6.307155e-06 0.0000693787
## 2 CE_16:0.1 0.9155465 2.225613e-16 4.029935 5.853538e-05 0.0003219446
## 3 TG_52:3 -0.8818559 -5.731354e-17 -3.881640 1.081696e-04 0.0003904006
## 4 TG_53:7 -0.8665732 5.344425e-17 -3.814371 1.419638e-04 0.0003904006
## 5 DG_36:2 -0.8525358 4.763723e-17 -3.752583 1.815657e-04 0.0003994444
## 6 TG_52:2.1 -0.8391009 7.714598e-17 -3.693447 2.290206e-04 0.0004198711
## 7 TG_53:2 -0.7702227 3.874216e-17 -3.390267 7.160239e-04 0.0011251804
## 8 DG_34:1 -0.7479340 -1.896877e-16 -3.292160 1.016891e-03 0.0013982253
## 9 TG_54:2 -0.7253484 -4.478293e-17 -3.192745 1.437954e-03 0.0017574991
## 10 Cholesterol 0.6866663 3.254021e-16 3.022480 2.548805e-03 0.0028036859
## 11 TG_53:6 -0.4913213 1.675669e-17 -2.162635 3.072386e-02 0.0307238648
## B
## 1 3.5407333
## 2 1.4405388
## 3 0.8669213
## 4 0.6138069
## 5 0.3852174
## 6 0.1699332
## 7 -0.8800780
## 8 -1.2006081
## 9 -1.5158089
## 10 -2.0331984
## 11 -4.2129728
visualize
full_results_final %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>%
ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
significant_final<-full_results_final[full_results_final$adj.P.Val<p_cutoff,]
head(significant_final)
## symbols logFC AveExpr t P.Value adj.P.Val
## 1 TG_53:3 -1.0295330 3.430077e-17 -4.531666 6.307155e-06 0.0000693787
## 2 CE_16:0.1 0.9155465 2.225613e-16 4.029935 5.853538e-05 0.0003219446
## 3 TG_52:3 -0.8818559 -5.731354e-17 -3.881640 1.081696e-04 0.0003904006
## 4 TG_53:7 -0.8665732 5.344425e-17 -3.814371 1.419638e-04 0.0003904006
## 5 DG_36:2 -0.8525358 4.763723e-17 -3.752583 1.815657e-04 0.0003994444
## 6 TG_52:2.1 -0.8391009 7.714598e-17 -3.693447 2.290206e-04 0.0004198711
## B
## 1 3.5407333
## 2 1.4405388
## 3 0.8669213
## 4 0.6138069
## 5 0.3852174
## 6 0.1699332
dim(significant_final)
## [1] 11 7
lets observe the -log10(P.value) of the significant lipids with a plot
graph<-ggplot(data=significant_final, aes(x=symbols, y=-log10(P.Value))) +
geom_bar(stat="identity", fill="steelblue")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
graph
observe the lipids logFC of the lipid
graph_FC<-ggplot(data=significant_final, aes(x=symbols, y=logFC)) +
geom_bar(stat="identity", fill="steelblue")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
graph_FC
perform wilcoxon test on the significant lipids
#extract the significant lipids
lipid_stat<-lipid_scaled[,colnames(lipid_scaled) %in% significant_final$symbols ]
number.of.samples<-ncol(lipid_stat)
# Create a output
p.value.vector <- c()
# Assign it to a list
p.value.list <- numeric(number.of.samples)
# Add the sample type and preprocess the dataframe
lipid_stat<-as.data.frame(cbind(lipid_stat,group=clinical_data$NAFLD))
lipid_steatosis0<-lipid_stat[lipid_stat$group=="steatosis0", ]
lipid_steatosis1<-lipid_stat[lipid_stat$group=="steatosis1", ]
lipid_steatosis0<-t(lipid_steatosis0)
lipid_steatosis1<-t(lipid_steatosis1)
input.combined.b <- cbind(lipid_steatosis0[,1:120],lipid_steatosis1[,1:21])
input.combined.b<-as.data.frame(input.combined.b)
input.combined.b<-input.combined.b[-12,]
input.combined.b = data.frame(lapply(input.combined.b, function(x) as.numeric(as.character(x))),check.names=F, row.names = rownames(input.combined.b))
# Establishing variable with number of observations
number.of.samples <-nrow(input.combined.b)
# Establishing vector with number of observations
vector.number.of.samples <-1:number.of.samples
for (i in vector.number.of.samples) {
steatosis0 <- unlist(input.combined.b[i, 1:120])
steatosis1 <- unlist(input.combined.b[i, 121:141])
wilcox.test.i <- wilcox.test(steatosis0, steatosis1,
mu=0,
alt="two.sided",
paired=F,
conf.int=F,
conf.level=0.95,
exact=F
)
p.value.list[i] <- wilcox.test.i["p.value"]
}
create a dataframe
#lipid list names to make a summary
lipid.list.df <- paste(rownames(input.combined.b) )
lipid.list.df<-as.data.frame(lipid.list.df)
p.value<- data.frame(matrix(unlist(p.value.list ), nrow=length(p.value.list ), byrow=TRUE))
colnames(p.value)<-"p.value"
#binding lipid list and pvalue to create summary list
wilcox.test.results <- cbind(lipid.list.df, p.value)
print(wilcox.test.results)
## lipid.list.df p.value
## 1 CE_16:0.1 7.555428e-05
## 2 Cholesterol 1.376790e-03
## 3 DG_34:1 5.172790e-04
## 4 DG_36:2 9.612430e-05
## 5 TG_52:2.1 5.364438e-05
## 6 TG_52:3 3.013631e-05
## 7 TG_53:2 1.296171e-04
## 8 TG_53:3 2.821378e-06
## 9 TG_53:6 3.473401e-03
## 10 TG_53:7 2.849056e-05
## 11 TG_54:2 8.942429e-05
print(wilcox.test.results$lipid.list.df)
## [1] "CE_16:0.1" "Cholesterol" "DG_34:1" "DG_36:2" "TG_52:2.1"
## [6] "TG_52:3" "TG_53:2" "TG_53:3" "TG_53:6" "TG_53:7"
## [11] "TG_54:2"
visualize the p values from wilxocon test
bar_graph<-ggplot(data=wilcox.test.results, aes(x=lipid.list.df, y=-log10(p.value))) +
geom_bar(stat="identity", fill="steelblue")+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
#theme_minimal()
bar_graph
Order the clinical data
clinical_data_boxplot<-clinical_data[match(colnames(input.combined.b), clinical_data$sample_ID),]
boxplot_data<-t(input.combined.b)
boxplot_data<-as.data.frame(boxplot_data)
boxplot_data<-add_column(boxplot_data, samples = clinical_data_boxplot$NAFLD, .before = 1)
head(boxplot_data)
## samples CE_16:0.1 Cholesterol DG_34:1 DG_36:2 TG_52:2.1 TG_52:3
## N1 steatosis0 1.16329901 0.18285982 -1.251192 -1.2132063 -0.8970818 -1.222416
## N2 steatosis0 0.43946354 -0.19793041 -1.246836 -0.8620915 -0.8583501 -1.315292
## N3 steatosis0 0.03958326 -0.30009357 -1.452855 -1.2032769 -0.9457813 -1.262850
## N4 steatosis0 0.16038153 -0.05679797 -1.336211 -0.8620326 -0.9107814 -1.199429
## N5 steatosis0 0.33983494 0.13088503 -1.318653 -0.8693971 -0.8836531 -1.241908
## N6 steatosis0 0.49945675 0.73867030 -1.389051 -1.3103400 -0.9648729 -1.208322
## TG_53:2 TG_53:3 TG_53:6 TG_53:7 TG_54:2
## N1 -1.2597033 -1.2412499 -1.251368 -1.057828 -0.9565674
## N2 -0.9259357 -0.5980203 -1.251368 -1.057828 -0.7283276
## N3 -1.2597033 -0.9276351 -1.251368 -1.057828 -0.9942757
## N4 -1.0089002 -1.2412499 -1.251368 -1.057828 -0.7656261
## N5 -0.9301686 -0.8299763 -1.251368 -1.004041 -0.7705726
## N6 -1.2597033 -0.9003478 -1.251368 -1.057828 -0.8283395
#write.csv(boxplot_data,"lipids_for_analysis.csv")
create a boxplot for visualizing the difference between the steatosis-0 and steatosis-1 samples.
boxplot_data$samples=as.factor(boxplot_data$samples)
my_lipid_list <- c("CE_16:0.1","Cholesterol","DG_34:1","DG_36:2","TG_52:2.1","TG_52:3","TG_53:2","TG_53:3","TG_53:6","TG_53:7","TG_54:2")
my_plot_list <- vector(mode = "list", length = 11)
for (i in 1:11) {
my_comparisons <- list( c("steatosis0", "steatosis1") )
p <-ggboxplot(boxplot_data, x = "samples", y = my_lipid_list[i],
color = "samples", palette = "jco",add = "jitter")+
stat_compare_means(label.y = 16)
my_plot_list[[i]] <- p
}
print(my_plot_list[1:11])
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
create complex heatmap with all the samples in a single plot
ht1 = Heatmap(input.combined.b[,1:120], name = "steatosis0", row_title = "", column_title = "steatosis0")
ht2 = Heatmap(input.combined.b[121:141], name = "steatosis1", row_title = "", column_title = "steatosis1")
ht_list = ht1 + ht2
draw(ht_list, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = "Comparison between steatosis0 and steatosis1", column_title_side = "bottom")
lets annotate the samples along with their gender to view their differences in heatmap
clinical_data_F<-subset(clinical_data, sex=="F",select=sample_ID:sex)
clinical_data_M<-subset(clinical_data, sex=="M",select=sample_ID:sex)
#lets add their gender
clinical_data_F$sample_sex<-sub("^","F.",clinical_data_F$sample_ID)
clinical_data_M$sample_sex<-sub("^","M.",clinical_data_M$sample_ID)
clinical_data_both_sex<-rbind(clinical_data_F,clinical_data_M)
#order the newly created metadata to the colnames of the lipid expression data
clinical_data_both_sex_ordered<-clinical_data_both_sex[match(colnames(input.combined.b), clinical_data_both_sex$sample_ID),]
lets check if the column names of the expression lipid data is similar to the clinical data
all(colnames(input.combined.b) == clinical_data_both_sex_ordered$sample_ID)
## [1] TRUE
head(input.combined.b[,1:10])
## N1 N2 N3 N4 N5 N6
## CE_16:0.1 1.1632990 0.4394635 0.03958326 0.16038153 0.3398349 0.4994568
## Cholesterol 0.1828598 -0.1979304 -0.30009357 -0.05679797 0.1308850 0.7386703
## DG_34:1 -1.2511917 -1.2468364 -1.45285479 -1.33621146 -1.3186531 -1.3890511
## DG_36:2 -1.2132063 -0.8620915 -1.20327690 -0.86203257 -0.8693971 -1.3103400
## TG_52:2.1 -0.8970818 -0.8583501 -0.94578129 -0.91078143 -0.8836531 -0.9648729
## TG_52:3 -1.2224155 -1.3152917 -1.26284959 -1.19942874 -1.2419081 -1.2083217
## N7 N8 N9 N10
## CE_16:0.1 -0.7594497 -1.0764997 -0.50037377 0.06703728
## Cholesterol -0.9287608 -0.8771608 -0.33348685 0.08801765
## DG_34:1 -1.1093491 -0.5819656 -0.80144507 -1.26994037
## DG_36:2 -1.1479058 0.3503971 -0.79397421 -0.97467791
## TG_52:2.1 -0.8771772 -0.6298855 -0.75281917 -0.92615116
## TG_52:3 -0.9137006 -0.7715755 -0.08312642 -1.22462884
now lets replace the column names of the expression data with the newly created sample_ID+gender column
colnames(input.combined.b)<-clinical_data_both_sex_ordered$sample_sex
head(input.combined.b[,1:10])
## F.N1 F.N2 F.N3 F.N4 F.N5 F.N6
## CE_16:0.1 1.1632990 0.4394635 0.03958326 0.16038153 0.3398349 0.4994568
## Cholesterol 0.1828598 -0.1979304 -0.30009357 -0.05679797 0.1308850 0.7386703
## DG_34:1 -1.2511917 -1.2468364 -1.45285479 -1.33621146 -1.3186531 -1.3890511
## DG_36:2 -1.2132063 -0.8620915 -1.20327690 -0.86203257 -0.8693971 -1.3103400
## TG_52:2.1 -0.8970818 -0.8583501 -0.94578129 -0.91078143 -0.8836531 -0.9648729
## TG_52:3 -1.2224155 -1.3152917 -1.26284959 -1.19942874 -1.2419081 -1.2083217
## F.N7 F.N8 F.N9 F.N10
## CE_16:0.1 -0.7594497 -1.0764997 -0.50037377 0.06703728
## Cholesterol -0.9287608 -0.8771608 -0.33348685 0.08801765
## DG_34:1 -1.1093491 -0.5819656 -0.80144507 -1.26994037
## DG_36:2 -1.1479058 0.3503971 -0.79397421 -0.97467791
## TG_52:2.1 -0.8771772 -0.6298855 -0.75281917 -0.92615116
## TG_52:3 -0.9137006 -0.7715755 -0.08312642 -1.22462884
extract the file used to create the heatmap
file_heatmap<-t(input.combined.b)
file_heatmap<-as.data.frame(file_heatmap)
all(rownames(file_heatmap) == clinical_data_both_sex_ordered$sample_sex)
## [1] TRUE
file_heatmap<-add_column(file_heatmap,sample=clinical_data_both_sex_ordered$NAFLD,.before=1)
#head(file_heatmap)
#dim(file_heatmap)
#write.csv(file_heatmap,"lipid_analysis.csv")
view the heatmap for significant lipids by subsetting the samples into smaller categories.
ht1 = Heatmap(input.combined.b[,1:30], name = "steatosis0", row_title = "", column_title = "steatosis0")
ht2 = Heatmap(input.combined.b[,31:60], name = "steatosis0", row_title = "", column_title = "steatosis0")
ht3 = Heatmap(input.combined.b[,61:90], name = "steatosis0", row_title = "", column_title = "steatosis0")
ht4 = Heatmap(input.combined.b[,91:120], name = "steatosis0", row_title = "", column_title = "steatosis0")
draw(ht1, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = " steatosis0", column_title_side = "bottom")
draw(ht2, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = " steatosis0 ", column_title_side = "bottom")
draw(ht3, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = " steatosis0", column_title_side = "bottom")
draw(ht4, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = " steatosis0", column_title_side = "bottom")
ht5 = Heatmap(input.combined.b[,121:141], name = "steatosis1", row_title = "", column_title = "steatosis1")
draw(ht5, row_title = "Lipids", row_title_gp = gpar(col = "black"),
column_title = " steatosis1", column_title_side = "bottom")
lets create a correlation network
set.seed(111)
Data2<-t(input.combined.b)
corMat <- cor(Data2) # Correlate data
head(corMat)
## CE_16:0.1 Cholesterol DG_34:1 DG_36:2 TG_52:2.1 TG_52:3
## CE_16:0.1 1.0000000 0.8500356 -0.7254860 -0.7240950 -0.7718605 -0.5241541
## Cholesterol 0.8500356 1.0000000 -0.5952215 -0.6494509 -0.6631440 -0.3506181
## DG_34:1 -0.7254860 -0.5952215 1.0000000 0.8576533 0.9495881 0.7917609
## DG_36:2 -0.7240950 -0.6494509 0.8576533 1.0000000 0.8757690 0.6876598
## TG_52:2.1 -0.7718605 -0.6631440 0.9495881 0.8757690 1.0000000 0.7471021
## TG_52:3 -0.5241541 -0.3506181 0.7917609 0.6876598 0.7471021 1.0000000
## TG_53:2 TG_53:3 TG_53:6 TG_53:7 TG_54:2
## CE_16:0.1 -0.6755248 -0.7026092 -0.5763073 -0.6052691 -0.6623299
## Cholesterol -0.5533712 -0.5618185 -0.4802842 -0.4557635 -0.5734958
## DG_34:1 0.8932375 0.8427915 0.8708539 0.8297891 0.8073444
## DG_36:2 0.7601000 0.7698806 0.7926969 0.7282854 0.8003767
## TG_52:2.1 0.9030812 0.8679832 0.8465218 0.8285776 0.8019836
## TG_52:3 0.7377178 0.7564427 0.8649190 0.9224849 0.6230012
plot the graph
Groups <- c("Cholesteryl ester","Cholesterol",rep("Diacylglycerol",2),rep("Triacylglycerol",7))
Graph_pcor <- qgraph(corMat, graph = "pcor", layout = "spring",theme = "colorblind", threshold = "BH",
sampleSize = nrow(Data2), alpha = 0.05,groups=Groups,
legend.cex = 0.55, vsize = 10,
esize = 15, borders = T,label.fill.vertical=0.3,
label.scale=T,edge.labels=T,edge.label.bg=T,repulsion=0.35)
lets obtain the CSV files of the significant data for the random forest model
#write.csv(clinical_data,'clinical_data_significant.csv')
#write.csv(input.combined.b,'lipid_data_significant.csv')
obtain the session info
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] igraph_1.2.6 qgraph_1.6.9
## [3] limma_3.46.0 ComplexHeatmap_2.6.2
## [5] ggpubr_0.4.0 pmp_1.2.1
## [7] EnhancedVolcano_1.8.0 ggrepel_0.9.1
## [9] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [11] Biobase_2.50.0 MatrixGenerics_1.2.1
## [13] matrixStats_0.59.0 GenomicRanges_1.42.0
## [15] GenomeInfoDb_1.26.7 IRanges_2.24.1
## [17] S4Vectors_0.28.1 BiocGenerics_0.36.1
## [19] rstatix_0.7.0 forcats_0.5.1
## [21] stringr_1.4.0 dplyr_1.0.6
## [23] purrr_0.3.4 readr_1.4.0
## [25] tidyr_1.1.3 tibble_3.1.1
## [27] ggplot2_3.3.3 tidyverse_1.3.1
## [29] readxl_1.3.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 tidyselect_1.1.1 RSQLite_2.2.7
## [4] AnnotationDbi_1.52.0 htmlwidgets_1.5.3 BiocParallel_1.24.1
## [7] munsell_0.5.0 codetools_0.2-16 withr_2.4.2
## [10] colorspace_2.0-1 highr_0.9 ggalt_0.4.0
## [13] knitr_1.33 rstudioapi_0.13 ggsignif_0.6.1
## [16] Rttf2pt1_1.3.8 GenomeInfoDbData_1.2.4 mnormt_2.0.2
## [19] bit64_4.0.5 vctrs_0.3.8 generics_0.1.0
## [22] xfun_0.23 itertools_0.1-3 randomForest_4.6-14
## [25] R6_2.5.0 ggbeeswarm_0.6.0 clue_0.3-59
## [28] locfit_1.5-9.4 bitops_1.0-7 cachem_1.0.5
## [31] DelayedArray_0.16.3 assertthat_0.2.1 scales_1.1.1
## [34] nnet_7.3-14 beeswarm_0.4.0 gtable_0.3.0
## [37] ash_1.0-15 Cairo_1.5-12.2 rlang_0.4.10
## [40] genefilter_1.72.1 GlobalOptions_0.1.2 splines_4.0.3
## [43] extrafontdb_1.0 impute_1.64.0 broom_0.7.6
## [46] checkmate_2.0.0 yaml_2.2.1 reshape2_1.4.4
## [49] abind_1.4-5 modelr_0.1.8 backports_1.2.1
## [52] Hmisc_4.5-0 extrafont_0.17 tools_4.0.3
## [55] lavaan_0.6-9 psych_2.1.6 ellipsis_0.3.2
## [58] jquerylib_0.1.4 RColorBrewer_1.1-2 Rcpp_1.0.6
## [61] plyr_1.8.6 base64enc_0.1-3 zlibbioc_1.36.0
## [64] RCurl_1.98-1.3 rpart_4.1-15 pbapply_1.4-3
## [67] GetoptLong_1.0.5 haven_2.4.1 cluster_2.1.0
## [70] fs_1.5.0 magrittr_2.0.1 data.table_1.14.0
## [73] openxlsx_4.2.3 circlize_0.4.13 reprex_2.0.0
## [76] pcaMethods_1.82.0 tmvnsim_1.0-2 hms_1.1.0
## [79] evaluate_0.14 xtable_1.8-4 XML_3.99-0.6
## [82] rio_0.5.26 jpeg_0.1-8.1 gridExtra_2.3
## [85] shape_1.4.6 compiler_4.0.3 maps_3.3.0
## [88] KernSmooth_2.23-17 crayon_1.4.1 htmltools_0.5.1.1
## [91] corpcor_1.6.9 Formula_1.2-4 geneplotter_1.68.0
## [94] lubridate_1.7.10 DBI_1.1.1 dbplyr_2.1.1
## [97] proj4_1.0-10.1 MASS_7.3-53 Matrix_1.2-18
## [100] car_3.0-10 cli_2.5.0 pkgconfig_2.0.3
## [103] foreign_0.8-80 xml2_1.3.2 foreach_1.5.1
## [106] pbivnorm_0.6.0 annotate_1.68.0 vipor_0.4.5
## [109] bslib_0.2.5.1 missForest_1.4 XVector_0.30.0
## [112] rvest_1.0.0 digest_0.6.27 rmarkdown_2.8
## [115] cellranger_1.1.0 htmlTable_2.2.1 curl_4.3.1
## [118] gtools_3.8.2 rjson_0.2.20 glasso_1.11
## [121] lifecycle_1.0.0 nlme_3.1-149 jsonlite_1.7.2
## [124] carData_3.0-4 fansi_0.4.2 pillar_1.6.1
## [127] lattice_0.20-41 ggrastr_0.2.3 fastmap_1.1.0
## [130] httr_1.4.2 survival_3.2-7 glue_1.4.2
## [133] zip_2.1.1 fdrtool_1.2.16 png_0.1-7
## [136] iterators_1.0.13 bit_4.0.4 stringi_1.5.3
## [139] sass_0.4.0 blob_1.2.1 latticeExtra_0.6-29
## [142] memoise_2.0.0